disp('hqlfc');

Nc = 40;

vh = linspace(1e-4, 1e-1, 80);
vres = zeros(1, length(vh));



h0 = 0.1;

vk = linspace(0, 10, 800);
mr = zeros(Nc, length(vk));

for ik = 1:length(vk)
    
    k = vk(ik);
    mr(:,ik) = fn_solve_ft(k, Nc, h0);
    
end

% close(figure(1));
figure(1);

for n = 1:Nc
    vr = mr(n, :);
    figure(1); hold on;
    plot(vk, vr, '-o');
end

mr2 = zeros(size(mr));

for n = 1:Nc
    
    vr = mr(n, :);
    vplus = exp(-n/2*abs(vk));
    vr2 = vr.*vplus;
    mr2(n, :) = vr2;
    
end

% close(figure(2));
figure(2);

for n = 1:Nc
    vr = mr2(n, :);
    figure(2); hold on;
    plot(vk, vr, '-o');
end

% close(figure(4));
figure(4);
tv = sum(mr2);
figure(4); hold on;
plot(vk, tv, '-o');


for ih = 1:length(vh)
    
    h0 = vh(ih);
    
    vk = linspace(0, 10, 800);
    mr = zeros(Nc, length(vk));
    
    for ik = 1:length(vk)
        
        k = vk(ik);
        mr(:,ik) = fn_solve_ft(k, Nc, h0);
        
    end
    
    
    % close(figure(1));
    % figure(1);
%     for n = 1:Nc
%         
%         figure(1); hold on;
%         plot(vk, mr(n,:), '-o');
%         
%     end
    
    mr2 = zeros(size(mr));
    
    for n = 1:Nc
        
        vr = mr(n, :);
        vplus = exp(-n/2*abs(vk));
        vr2 = vr.*vplus;
        mr2(n, :) = vr2;
        
    end
    
%     for n = 1:Nc
%         
%         figure(2); hold on;
%         plot(vk, mr2(n,:), '-o');
%         
%     end
    
    vr2 = sum(mr2);
    
    figure(3); hold on; 
    plot(vk, vr2, '-o');
    
    vk_nt = zeros(size(vk));
    vlen = vk(2:end) - vk(1:end-1);
    vk_nt(1:end-1) = 1/2*vlen;
    vk_nt(2:end) = vk_nt(2:end) + 1/2*vlen;
    
    disp(sum(vk_nt.*vr2));
    
    vres(ih) = sum(vk_nt.*vr2);
    
end

figure(999); hold on;
plot(vh, vres, '-o');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function res = fn_solve_ft(k, Nc, h0)

Vn = fn_Vn(k, Nc);
An = fn_An(h0, Nc);
T = fn_T(k, Nc);

M = eye(Nc) + T*diag(An);

fprintf('cond(M) = %f \n', cond(M));

res = M\Vn;

end


function Vn = fn_Vn(k, Nc)

Vn = zeros(Nc, 1);

for n = 1:Nc
    
    Vn(n) = exp(-n/2*abs(k));
    
end

end


function An = fn_An(h0, Nc)

An = zeros(1, Nc);

for n = 1:Nc
    
    An(n) = sinh(h0/2)^2/sinh(h0/2*(n+1))^2;
    
end

end


function T = fn_T(k, Nc)

T = zeros(Nc, Nc);

if abs(k)<1e-10
    
    for n = 1:Nc
        for m = 1:Nc
            T(n,m) = 2*min(n,m);
        end
    end
    
else
    
    for n = 1:Nc
        for m = 1:Nc
            T(n,m) = (1+exp(-abs(k)))/(1-exp(-abs(k))) ...
                *(exp(-abs(n-m)/2*abs(k)) - exp(-abs(n+m)/2*abs(k)));
        end
    end
    
end

end